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ABSTRACT 


Daily  observations  of  sea-surface  temperature  (SST)  taken  just  north 
of  Pt.  Sur,  California  over  the  past  12  years  are  examined  for  statistical  and 
oceanographic  content.  Spectral  analysis  of  the  data  reveals  an  annual  cycle 
and  two  harmonics,  but  no  obvious  influence  from  the  local  winds  or  tides.  A 
statistical  model  has  been  constructed  based  on  a  decomposition  of  the  data 
into  its  various  deterministic  and  random  components.  As  expected,  the  model 
provides  forecasts  which  are  significantly  better  than  those  that  would  be 
obtained  from  a  simple  climatological  treatment  of  the  data  alone.  The  model 
is  also  used  to  generate  simulated  series,  which  in  turn  are  used  as  a  basis 
to  estimate  adequate  sampling  rates. 

In  6  years  of  the  12  year  record,  an  abrupt  major  decrease  in  SST  was 
observed  between  February  and  April  and  has  been  identified  as  the  spring 
transition  to  coastal  upwelling.  One  such  transition  was  found  to  occur 
almost  simultaneously  at  three  locations,  including  Pt.  Sur,  along  the 
California  coast  between  Pt.  Conception  and  San  Francisco.  The  data  also 
indicate  that  a  well-defined  annual  temperature  cycle  does  exist  at  Pt.  Sur 
although  Its  range  is  reduced  and  its  maximum  occurs  later  in  the  year  than 
would  normally  be  expected  at  this  latitude.  Three  El  Nino  events  are  clearly 
evident  In  the  12-year  record  including  the  present  El  Nino;  SSTs  have 
increased  by  almost  3C  at  Pt.Sur  as  a  result  of  the  present  El  Nino  warming. 
Over  the  whole  12  year  period  an  average  increase  in  temperature  of  1.6C  was 
found. 
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ANALYSIS  OF  A  12-YEAR  RECORD  OF  SEA-SURFACE 
TEMPERATURES  OFF  PT.  SUR,  CALIFORNIA 


I.  Introduction 

Sea-surface  temperatures  (SSTs)  are  collected  on  a  regular  basis  (SIO 
Reference  81-30)  along  the  California  coast  at  approximately  2b  locations.  In 
view  of  the  increasing  interest  in  coastal  circulation  with  respect  to  the 
dispersal  of  pollutants  and  a  continuing  interest  in  fisheries  biology,  it  is 
surprising  that  few  systematic  studies  have  been  conducted  to  examine  these 
records,  particularly  within  the  past  15  years.  Perhaps  the  availability  of 
high-quality  tenperature  data  along  the  continental  shelves  from  new  and 
sophisticated  instruments  has  caused  us  to  overlook  this  resource  in  recent 
years.  However,  the  coastal  observations,  unlike  those  acquired  from  moored 
arrays,  often  extend  over  many  years  and  thus  provide  a  unique  opportunity  to 
examine  coastal  variability  over  relatively  long  periods.  For  locations  where 
the  measuring  site  has  a  good  exposure  to  the  adjacent  continental  shelf  and 
slope,  measurements  of  SST  may  be  particularly  revealing  with  respect  to  some 
of  the  physical  processes  that  occur  regionally  as  well  as  locally.  In  this 
regard  a  12-year  record  of  SSTs  off  the  central  California  coast  at  Granite 
Canyon ,  15km  south  of  Monterey,  has  been  chosen  for  the  initial  study.  These 
data  have  been  acquired  daily  by  the  Marine  Culture  Laboratory  of  the 
California  Fish  and  Game  Commission  since  March  of  1971.  The  Granite  Canyon 
site  where  the  observations  are  made  (see  Fig.  1)  has  an  excellent  exposure  to 
the  deep  ocean  with  the  continental  shelf  extending  less  than  10km  offshore. 

Previous  investigations  along  the  U.S.  West  Coast  Indicate  that  the 
oceanographic  information  content  of  coastal  SST's  is  usually  high. 

SSTs  have  been  shown  (a)  to  be  useful  indicators  of  ocean  temperature 
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FIGURE  1 


Measurement  Locations.  The  principal  time  series  considered  here 
was  taken  at  Granite  Canyon.  Qualitative  comparisons  are  made  to 
fragmentary  series  of  subsurface  temperatures  close  by,  and  to  a 
short  time  series  at  Diablo  Canyon.  Comparisons  of  the  Granite 
Canyon  SSTs  will  be  made  to  extensive  time  series  of  SSTs  at  the 
Hopkins  Marine  Station  at  Pacific  Grove  and  to  SSTs  at  the 
Farallon  Islands. 


variability,  (b)  to  be  often  representative  of  phenomena  occurring  over  wide 
regions,  (c)  to  be  related  to  other  ocean  and  atmospheric  variables,  and  (d) 
to  have  consistent  internal  structures.  Based  on  surface  and  subsurface 
temperatures  along  the  California  coast,  Reid,  Roden  and  Wyllie  (1958) 
indicated  two  different  causes  of  seasonal  variability  in  SSTs.  The  first  is 
coastal  upwelling.  North  of  Pt.  Conception  the  upwelling  process  decreases 
the  seasonal  range  in  temperature  and  increases  the  duration  of  cooler 
temperatures.  The  second  cause  stems  from  the  importance  of  the  coastal 
countercurrent.  This  flow  to  the  northwest  transports  warm  water  from  the 
south  along  the  California  coast,  increasing  surface  temperatures  during  the 
fall.  Also,  Reid,  Roden  and  Wyllie  (1958)  found  that  variations  in 
temperature  along  the  U.S.  West  Coast  show  considerable  coherence  over  great 
distances  and  that  major  anomalies  in  temperature  persist  for  several  months. 

Using  coastal  SSTs  (and  sea  level)  Stewart,  Zetler,  and  Taylor  (1958) 
verified  the  existence  of  anomalous  conditions  off  the  U.S.  West  Coast  during 
1957  and  1958.  Robinson  (1957)  Indicated  the  value  and  representativeness  of 
coastal  SSTs.  She  found  that  coastal  SSTs  and  air  temperatures  were  highly 
correlated  and  that  SST  anomalies  identified  at  coastal  locations  could  be 
correlated  with  SST  anomalies  up  to  250km  offshore.  Roden  (1963)  found  that 
non-annual  extreme  temperature  fluctuations  were  coherent  and  in  phase  over 
distances  of  200-300km.  He  also  found  that  the  average  persistence  of 
non-annual  minimum  temperatures  was  3-4  months  while  the  average  persistence 
of  non-annual  maximum  tenperatures  was  3-5  months  along  the  coast  of  northern 
California.  Spectral  analysis  Indicated  the  presence  of  an  annual  cycle  only, 
however. 

More  recently  List  and  Koh  (1976)  analyzed  SSTs  along  the  California 
coast  over  several  years  using  spectral  analysis  and  found  that  many  of  the 
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coastal  stations  exhibited  strong  spectral  peaks  at  14.7  days  and  that  these 
peaks  were  probably  the  result  of  tidal  aliasing  and/or  the  spring  and  neap 
tides.  They  also  found  that  a  well-defined  seasonal  variation  in  SST  was 
essentially  absent  between  Pt.  Conception  and  Monterey  Bay  and  that  certain 
fluctuations  occurred  simultaneously  at  most  coastal  stations. 

Surface  and  subsurface  tenperatures  across  Monterey  Bay  were  collected  on 
a  regular  basis  by  the  Hopkins  Marine  Station  from  1929  to  1974.  Results  of 
these  measurements  have  been  reported  by  Skogsberg  (1936),  Skogsberg  and 
Phelps  (1946),  and  Bolin  and  Abbott  (1963).  Much  of  the  data  presented  in 
these  references  were  time-averaged  providing  monthly  mean  thermal  conditions 
over  the  Bay.  Skogsberg  used  these  data  to  define  three  separate  oceanic 
seasons  for  the  Monterey  Bay  area. 

Our  purpose  in  this  paper  is  threefold:  first,  to  model  (i.e.  describe 
mathematically  and  statistically)  the  behavior  of  the  12  year  Granite  Canyon 
data  set;  second,  to  use  the  model  and  other  statistical  techniques  to  project 
or  predict  the  data  and  estimate  time  scales  of  variability  in  SSTs,  and 
third,  to  provide  a  descriptive  interpretation  of  the  Granite  Canyon  data  from 
an  oceanographic  viewpoint. 

2.  The  Granite  Canyon  Data 

SSTs  at  Granite  Canyon  are  measured  dally  at  approximately  0800 
local  time.  The  raw  data  are  shown  in  Figure  2.  Because  the  data  are  taken 
only  once  every  24  hours,  the  possibility  of  aliasing  by  the  predominant 
semidiurnal  tide  exists.  This  problem  will  be  addressed  later  in  the 
analysis.  As  mentioned  before  the  location  has  an  excellent  exposure  to  the 
deep  ocean.  Temperatures  are  read  to  the  nearest  0.1C  using  a  calibrated  VWR 
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RAW  SEA- SURFACE  TEMPERATURE  AT  GRANITE  CANYON-3/1 /71  TO  3/1/83 

MEASUREMENTS  TAKEN  AT  0800  PST 


FIGURE 


thermometer,  and  are  considered  accurate  to  about  ±  0.2C  (SIO  Reference  81-30). 
The  observations  span  12  years  starting  March  1,  1971  and  ending  on  February 
28,  1983.  SSTs  over  the  same  time  period  were  acquired  for  the  Hopkins  Marine 
Station  at  Pacific  Grove  and  the  Farallon  Islands  (see  Figure  1  for  location). 
The  accuracy  of  these  data  is  also  about  ±  0.2C. 

Since  the  data  from  Granite  Canyon  (GC)  have  only  been  available  since 
1971  and  since  they  have  not  been  compared  with  other  coastal  temperatures,  it 
was  of  interest  initially  to  consider  the  representativeness  of  these 
observations.  Consequently  we  have  compared  the  Granite  Canyon  data  with 
several  other  data  sets  for  periods  where  comparisons  could  be  made.  The  annual 
variations  in  SST  were  generally  similar  at  Granite  Canyon,  Diablo  Canyon  and 
the  Farallons,  but  temperatures  were  consistently  higher  at  Pacific  Grove 
during  the  spring  and  summer  months  than  they  were  at  the  other  three 
locations. 

Temperature  data  at  four  depths  between  113  and  510m  were  acquired  from  a 
moored  array  approximately  50km  south  of  Granite  Canyon  for  the  period  03/05/80 
to  04/07/80  and  compared  with  the  SST  data  at  Granite  Canyon,  see  Figure  3. 

The  period  was  chosen  because  it  includes  a  major  event  tentatively  identified 
as  the  spring  transition  to  upwelling,  a  phenomenon  which  has  been  recently 
identified  further  north,  along  the  coast  of  Oregon,  by  Huyer  et  al,  (1979). 
This  event  can  be  readily  identified  in  a  number  of  the  annual  cycles  in  the 
raw  data  (Figure  2,  and  Table  1).  The  drop  in  temperature  associated  with  this 
event  during  March  of  1980  can  also  be  seen  at  each  depth  in  Figure  3,  down  to 
510m. 

Figures  4,  5,  6  show  comparisons  of  SSTs  between  Granite  Canyon  and 
Diablo  Canyon,  between  Granite  Canyon  and  Pacific  Grove,  and  between  Granite 
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GRANITE  CANYON  AND  DIABLO  CANYON  COMPARED 


FIGURE  4 


GRANITE  CANYON  AND  THE  FARALLONS  COMPARED 


FIGURE 
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Canyon  and  the  Farallon  Islands  respectively  for  the  year  1980.  Diablo  Canyon 
Is  located  about  200  Km  south  of  Granite  Canyon  and  the  data  from  that  location 
were  recorded  at  5m  below  the  surface  on  the  continental  shelf,  just  off  the 
PG4E  nuclear  power  plant  located  on  the  adjacent  coast.  The  temperatures  at 
Pacific  Grove  are  collected  near  the  shore  in  a  sheltered  area  inside  Monterey 
Bay.  The  Farallon  Islands  data  contained  numerous  multi -day  gaps  which  were 
filled  by  linear  Interpolation.  The  same  spring  transition  shown  in  Figure  3 
is  seen  to  be  present  in  these  figures.  The  overall  patterns  in  the  annual 
cycle  at  these  four  locations  are  seen  to  be  generally  similar  at  least  during 
1980. 


3.  Gross  trends  in  the  Granite  Canyon  SST's 

We  discuss  here  gross  features  of  the  Granite  Canyon  SST  record.  These 
features  are  of  interest  per  se  and  are  also  a  prelude  to  the  more 
quantitative  decomposition  and  time  series  analysis  performed  in  the 
subsequent  sections. 

Figure  7  shows  a  scatter  plot  of  the  temperature  data  against  day  number. 
In  most  years  fewer  than  three  observations  were  missing  from  the  raw  data. 
Where  missing  data  were  found,  linearly  interpolated  values  were  substituted. 
Since  the  points  In  Figure  7,  unlike  those  in  Figure  3,  are  not  connected,  the 
quantization  of  the  data  (0.1C)  can  be  clearly  seen. 

The  most  interesting  feature  of  Figure  7  is  a  definite  upward  average 
tilt  to  the  data.  A  least-squares  fit  of  the  linear  trend  equation 
y(t)  *  o  +  it  to  the  data  gives  the  fit  y(t)=10.9  +  Q.00037t,  where  t  is  in 
days,  and  y(t)  Is  water  temperature  in  °C  on  day  t.  This  indicates  an 
"average"  temperature  rising  from  10.9  to  12.52C  over  the  twelve  year  measure¬ 
ment  period  (t  *  4380).  Note  that  this  linear  trend  is  simply  a  convenient  way 


11 


SEA-SURFACE 


(0  •*  CN  O 


o«  ni  3aniva3d^3i 


12 


9 


of  representing  the  overall  Increase  In  temperature  over  the  12-year  period. 
Higher  order  polynomial  fits  provide  only  slight  Improvement.  We  may,  how¬ 
ever,  be  observing  only  a  part  of  a  much  longer  cycle  of  fixed  or  variable 
period  which  would  be  revealed  in  a  much  longer  data  series.  The  statistical 

A 

significance  of  the  trend  parameter  6  *0.00037  (l.e.  whether  it  Is  really 
different  from  zero)  Is  complicated  by  the  presence  of  cycles  and  possible 
short  term  correlation  In  the  data.  Thus  standard  tests  that  6*0  based  on 
Independent  data  are  Invalid.  However,  using  simulation  techniques  described 
later,  we  have  found  the  standard  deviations  of  6  to  be  0.000011  and  hence 
reject  the  hypothesis  that  6  »  0  at  the  95%  level  of  confidence. 

The  gross  features  of  the  Granite  Canyon  SST's  clearly  Include  annual 
cycles  and  It  Is  apparent  from  Figures  2  and  7  that  these  annual  cycles 
contain  inter-annual  variability  as  well.  In  particular  the  oceanic  El  Ninos 
of  1972,  1976,  and  1983  are  clearly  vlsable.  The  possibility  of  a  temperature 
anomaly  In  1979  Is  also  vlsable.  Although  the  effects  of  El  Nino  events  are 
often  restricted  to  equatorial  regions  and  to  regions  along  the  South  American 
coast,  their  Influence  does  upon  occasion  extend  at  least  as  far  north  as 
central  California  (Bretschnelder  and  McLain,  1983). 

Underlying  patterns  In  the  data  can  be  exposed  by  performing  local 
smoothing.  (Since  the  computations  In  this  paper  were  done  using 
the  APL  language  (see  e.g.  Gilman  i  Rose,  1976)  with  an  experimental  APL 
Graphics  package,  GRAFSTAT3,  the  smoothing  was  performed  using  the  program 
FILTER  from  Ans combe  (1981,  p392).)  The  smoothing  computes  a  cosine-weighted 
moving  average  of  extent  m  to  give  a  smoothed  time  series  y(t)  : 
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m-1 

y(t)  »  I  y(t)w(t-j)  ,  (3.1) 

j=-(m-l) 

where  the  w(k)=  (l-cos[2ir(k+l)/(m+l)]),  for  k30,±l,±2,..  ,±(m-l),  (3.2) 

are  equal-spaced  ordinates  of  l-cos(2t).  The  gain  of  this  filter  is  given  by 

G(<d)~  - sin  U(m*l)2] -  . 

i»i(m+l)[l-<u2(m+l)2(2w)2] 

Details  are  given  in  Anscombe  (1981,  p.156). 

Smoothings  of  the  Granite  Canyon  SST's  for  various  values  of  (2m+l) 
together  with  the  number  of  filter  weights  used  are  given  in  Figure  8.  It  is 
clear  that  61  weights  (or  81  weights)  resolve  the  gross  details  adequately. 

It  is  interesting  to  note  that  there  is  not  only  a  sharp  drop  in  temperature 
In  the  spring,  but  also  a  rather  clearly  defined  and  relatively  abrupt  rise  in 
temperature  In  the  fall.  We  note  that  both  the  spring  and  fall  transitions 
occur  over  periods  relatively  short  in  comparison  with  the  annual  cycle. 

Since  well-defined  spring  transitions  could  be  clearly  identified  in  at 
least  6  out  of  the  12  years.  It  was  decided  to  tabulate  these  events  with 
respect  to  their  dates  of  occurrence,  duration,  and  associated  changes  in 
temperature.  The  results  are  shown  in  Table  1  below.  Because  the  spring 
transitions  add  complexity  to  the  annual  temperature  cycle,  the  harmonic 
content  of  this  data  is  undoubtedly  richer  by  their  presence.  We  will  say 
more  about  the  annual  variation  in  temperature  at  Granite  Canyon  in  Section 
10. 
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Table!* 


YEAR 

STARTING  DATE 

DURATION 

DELTA  T 

dT/dt 

1973 

March  9 

7  days 

13.0  -  10.0C 

0.43°/Day 

1974 

March  31 

5  days 

12.3  -  10.0 

0.46 

1977 

Feb.  22 

10  days 

13.1  -  9.3 

0.38 

1978 

April  25 

9  days 

13.0  -  8.5 

0.50 

1980 

March  11 

7  days 

14.5  -  10.4 

0.59 

1981 

March  23 

8  days 

13.1  -  10.0 

0.39 

*  A  lesser  transition  was  observed  on  December  19,  of  1981.  Whether  or  not 
this  event  corresponds  to  a  "spring"  transition  is  not  known. 


A  detailed  statistical  analysis  of  the  time  series  of  SSTs  will  now  be 
given.  This  will  be  used  in  Sections  8  and  9  to  give  a  more  detailed  exam¬ 
ination,  statistical  and  oceanographic,  of  the  SSTs. 

4.  Cyclic  fit  and  decomposition  of  the  data:  General  ideas 

We  proceed  to  decorqpose  the  SST  time  series  into  an  additive  model  of  the 

form 

y (t )  *  m(t)  +  s (t )  +  e(t ),  (4.1) 

where  the  trend  m(t)  can  be  defined  loosely  as  a  long  term  change  in  the  mean 
(Chatfleld,  1980,  p.13),  s (t )  consists  of  seasonal  and  cyclic  changes,  and 
e(t)  is  a  stationary  random  sequence  which  describes  the  irregular  fluctua¬ 
tions.  The  c(t)  sequence  is  assumed  to  have  mean  zero,  (i.e.  E(e(t))=  0), 
constant  variance,  and  possibly  some  serial  dependence.  That  m(t),  s(t)  and 
e(t)  are  treated  as  additive  implies  the  gross  assumption  that  this  dependency 
structure  in  the  e(t),  as  well  as  the  variance,  is  not  affected  by  m(t),  or 


The  model  (4.1)  is  essentially  a  conceptual  tool  for  thinking  with  and 
its  usefulness  and  approximate  validity  must  be  established  from  the  data. 


Thus  to  define  "scales"  in  the  data  these  different  components  must  be 
identified. 

The  term  m(t)  was  discussed  in  Section  3  and  here  we  let 

y x (t )  *  y (t )  -  m(t)  =  s(t)  +  e(t)  (4.2) 

insofar  as  it  can  be  observed.  With  the  estimate  m(t)  =  10.9  +  0.000374t 
removed  from  y(t)  we  denote  the  resultant  as  the  first  residual,  r^t), 
ignoring  for  the  present  possible  interaction  between  the  estimation  of  m(t), 
s(t)  and  e (t ). 

In  this  section  we  estimate  the  structure  of  s(t)  and  e(t)  by  analyzing 
rj(t)  .  Note  that  s(t)  includes  definite  seasonal  terms,  such 
as  the  obvious  annual  effect  in  E (y (t ) )  (and  its  harmonics),  and  terms  with 
periods  longer  than  a  year  which  might  be  interpreted  as  global  fits  to  the 
obvious  long-term  modulation  of  the  annual  cycle  which  can  be  seen  in  figures 
2  and  6.  Again  there  may  be  cycles  at  shorter  periods  associated  with  the 
tides  and  winds;  these  are  generally  difficult  to  distinguish  from 
quasi-cycles  "generated"  by  highly  correlated  errors  e(t).  Thus  we  assume 
that 


$(t) 


{yj  cos(uijt)  +  sin(u>jt)} 


(4.3) 


J 

■  I  A,  cos (w.t  +  9.)  ,  (4.4) 

j-1  J  J  J 
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where  Uj*  2»fj  *  2*/t\  are  (unknown)  frequencies  corresponding  to  periods  Tj, 

A.  *  (y?  +  is  the  amplitude,  e.»  tan~*(-e ./y.)  is  the  phase,  and 

J  Is  the  number  of  cycles  which  comprise  s(t).  Again  we  note  that  cycles 
with  periods  longer  than  1  year  may  represent  long-term  random  fluctuations. 

The  representation  (4.3)  of  the  component  s(t)  of  E(y(t))  as  a  sum  of 
deterministic  cycles  with  fixed  frequencies  is  a  convenience  only.  In 
particular  it  might  be  hazardous  to  predict  very  long  term  behavior  on  the 
basis  of  this  model  for  s(t)  and  m(t ).  Also,  while  our  model  may  identify 
cycles  with  periods  shorter  than  one  year,  such  cycles  may  not  represent  real 
cyclic  behavior,  but  rather  they  may  serve  only  to  accommodate  the 
non-sinusoidal  nature  of  the  year  cycle  and  the  very  sharp  transition  in  that 
cycle  that  often  occurs  during  spring  (see  Table  1).  Because  our  data  extends 
for  only  12  years,  it  will  be  impossible  to  detect  any  cycles  with  periods 
longer  than  half  that  length,  or  6  years. 

To  find  "hidden  periodicities"  in  s(t)  is  a  formidable  task,  even  if  the 
spectrum  of  the  random  component  e(t),  namely 

f(u>)  *  —  {1  +  2  l  p(k)  cos(ku)}  ,  0  <  u  <  u  ,  (4.S) 

"  k=l  “  “ 

where  p(j)  s  corr(e(t+j )e(t)),  were  constant  (flat)  i.e.  the  e(t)  are  "white 
noise".  (This  is  the  classical  problem  of  "hidden  periodicities"  going  back  to 
the  original  work  of  Schuster  (1898)  and  Fisher  (1929).)  To  find  these 
periodicities  in  the  presence  of  a  non-flat  spectrum  for  the  e(t)'s,  we  proceed 
as  follows: 

(1)  Examine  the  periodogram  of  the  first  residuals,  r^t),  to  find 
possibly  significant  periodicities.  Estimate  s(t)  as  s(t). 

Given  the  chosen  frequencies  u>j,  this  involves  obtaining  least 


squares  estimates  of  the  y  js  and  0j's  in  the  linear  model  (4.3). 
The  least  squares  estimates  of  the  amplitudes  A.  are  directly  pro- 

J 

portlonal  to  the  periodogram  values  at  the  given  frequencies. 

(2)  Examine  the  second  residuals  r2(t)=y (t)-s(t)-m(t)=r1(t)-s(t) 

to  determine  the  properties  of  e(t).  If  e(t)  can  re  represented, 
say,  as  a  pth-order  autoregressive  process, 

e(t)  a  aje(t-l)  +  ...  +  ape(t-p)  +  n(t),  (4.6) 

where  the  n(t)  are  uncorrelated.  Then  we  estimate  aj  through  ap 
by  least  squares  to  be  aj  through  ap  and  "pre-whiten"  the  first 
residuals,  r^t),  as 

r3(t)  ■  rj(t)  -  a^ft-l)  ...  -Spr^t-p).  (4.7) 

We  then  examine  the  periodogram  of  these  third  residuals  to 
ascertain  whether  the  s (t )  chosen  in  step  (1)  is  still  a 
"reasonable"  representation  of  the  seasonal  and  cyclic  changes 
(global  fit)  for  the  data.  The  "prewhitening"  is  performed  because 
highly  autocorrelated  data  will  elevate  spectral  values  near  zero 
frequency,  making  it  difficult  to  decide  whether  a  given  spectral 
peak  is  due  to  the  correlation  or  due  to  a  cycle  in  the  data.  Once 
the  autocorrelated  component  is  removed,  we  obtain  a  clearer  picture 

A 

of  what  cycles  should  actually  be  included  in  s (t ) . 

(3)  Using  r3(t)  to  point  out  a  possibly  refined  s(t),  re-examine  the 
structure  of  e(t)  via  the  (possibly  refined)  second  residuals, 

r2(t)^y(t)-m(t)-s(t).  Knowledge  of  the  structure  of  e(t)  is 
important  for  determing  "local  scales"  and  obtaining  the  "whitened" 
residuals  r^(t).  Thus,  if  it  is  determined  that  e(t)  has  pth 


order  autoregressive  structure,  possibly  different  from  that  in  (2) 
above,  confute  fourth  residuals 

r4(t )  =  r2(t)  -  ajr^t-l)  ...  -apr2(t-p).  (4.8) 

These  residuals  are  used  firstly  to  establish  the  overall  goodness 
of  fit  of  the  estimated  model,  since  the  r^(t)  process  will  have 
a  flat  spectrum  and  a  zero  level  correlogram,  (i.e. ,  "white  noise") 
if  the  model  is  reasonable.  Secondly  they  can  be  used  (e.g.  Tsiao  and 
Box,  1982)  to  establish  cross-correlations  of  the  Granite  Canyon  data 
with  other  SST  data.  Thirdly,  they  provide  an  estimated  distribution 
of  the  n(t)  which  will  be  useful  later  for  modelling  and  for  checking 
assumptions  during  statistical  testing. 

5.  Global  fit  and  cyclic  components 

Figure  11  shows,  in  the  top  panel,  a  cyclic  fit  s(t)  with  J  =  13  terms 

a 

to  the  first  residuals  (original  data  minus  long  term  linear  trend  m(t)).  The 
data,  since  leap  days  have  been  removed,  consists  of  N  -  4380  daily  observations 
of  temperature.  Thus  if  angular  frequency  is  re-expressed  as  discrete  frequencies 
of  the  form  wp  =  2-np/N,  where  p  is  the  "number  of  cycles  per  period  of  4380 
days",  then  the  cycles  removed  from  the  data  as  components  of  s(t),  have 
frequencies  p  *  12,24,94,  corresponding  to  periods  of  1  year,  6  months  and 
approximately  1.5  months,  with  the  remaining  p  =  2,3,4,5,6,7,8,9,11,18. 

The  determination  of  these  frequencies  follows  the  steps  in  Section  4. 

First,  Figure  10  shows  the  estimated  serial  correlations  of  the  first 
residuals 

N“K 

p  (k)  *  l  (y(t)  -  m(t))  (y(t+k)-m(t))/{(N-k)S2}  (5.1) 

t=l 

N— k 

*  l  r,(t)r,(t+k)  |(N-k)S2} 

t«l  1 
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FREQUENCY  IN  CYCLES 


TERM  CYCLIC  FIT  OF  LINEARLY  DETRENDED  DATA 


FIGURE  11 


1 


where  S2  Is  the  sample  variance  of  the  r^tj's.  This  correlogram  is 
dominated  by  the  yearly  cycle  and  therefore  can  not  be  used  to  postulate  a 
model  for  the  e(t),  or  to  Identify  more  than  the  yearly  cycle  itself. 

A  harmonic  spectral  analysis  of  the  p (k ) *s  is  given  by  the  periodogram 
of  the  first  residuals,  defined  as 

In(“p)  “  <An  (V  +  Bn  (wp^  /{2,rN)’  p  =  1-*’N/2  »  (5-2) 

where 

N  N 

AN(u»p)*  l  rx(t)  cos(tup),  and  BN(wp)  =  l  rx(t)  sin(tu>p), 

for  p  *  lf...»N/2#  are  the  real  and  complex  parts  of  the  discrete  Fourier 
transform  of  rj(t),  t=l,  ...  N.  Here  up  =  -j^-  .  Since  the  expected  value  of 

2 

the  periodogram  Is  (asymptotically),  for  uncorrelated  data,  o  /2n,  we  use 

a  normalized  periodogram.  Thus  the  normalized  periodogram  is  defined  to  be 

o 

the  result  (5.2)  multiplied  by  2n/S2,  where  S  is  the  sample  variance  of  the 
data  being  considered,  namely  r^t). 

The  logarithm  of  the  normalized  periodogram  of  the  first  residuals  from 
the  Granite  Canyon  data  Is  shown  in  Figure  9  with  confidence  bands  for 
individual  components  and  for  the  maximum  value  of  the  periodogram  components. 
These  are  based  on  the  approximation  that  each  value  in  the  normalized  period¬ 
ogram  is  independently  distributed  as  a  mean  one  exponential  random  variable 
(chi-square  with  2  degrees  of  freedom).  The  dominant  frequency  is  p  *  12, 
i.e. ,  a  period  of  one  year.  There  is,  in  the  expanded  scale  periodogram  in  the 
top  panel  of  Figure  9,  a  significant  harmonic  at  p  *  24  and  at  p  =  94.  This 
latter  frequency  Is  surprising,  since  it  is  not  exactly  a  harmonic  of  p  =  12. 
The  full  periodogram  in  the  lower  panel  suggests,  even  though  it  is  unsmoothed, 
an  autoregressive  type  spectrum  for  the  e(t).  Thus,  determination  of  cycles  at 
frequencies  below  12  Is  difficult.  Proceeding  to  the  prewhitening,  frequencies 
p  *  2,3,4,5,6,7,8,9,11  and  18  are  selected  for  removal,  in  addition  to  p  *  12, 
24  and  94. 


1 


24 


Removal  of  these  frequencies  gives  the  second  residuals  whose  periodogram 
is  shown  in  Figure  12.  On  the  expanded  scale,  note  the  zeros  where  cycles  have 
been  removed;  this  is  because  the  least  squares  estimates  of  the  amplitudes  in 
(4.4)  are  exactly  proportional  to  the  corresponding  periodogram  values.  Again 
note  in  the  lower  panel  of  Figure  12  the  autoregressive-like  spectrum.  (A 
spectrum  for  a  theoretical  AR(2)  process  using  parameters  computed  below  is 
shown  for  comparison.) 

Figure  13  shows  the  estimated  autocorrelations  of  the  second  residuals  as 
well  as  the  first  50  correlations  for  four  theoretical  AR(P)  models  computed 
from  the  first  four  estimated  correlations,  via  the  Yule-Walker  equations 
(Chatfield,  1980,  p.48).  Note  in  Figure  13  that  the  AR(P)models  of  order 
greater  than  2  do  not  give  a  much  better  fit  than  the  AR(2)  model  to  the 
estimated  autocorrelations. 

Notice  also  that  the  persistent  hump  at  about  21  days.  As  yet  we  have  no 
explanation  for  this  effect.  The  autocorrelations  in  Figure  12  are  extended  to 
show  lags  beyond  50  in  Figure  13,  but  there  is  no  indication  of  a  cycle  of 
period  21.  Such  a  cycle  would  show  up  as  a  wave  in  the  correlogram  with  period 
21,  not  simply  as  a  peak  at  lag  21  days. 

Table  3 

Estimated  correlation  coefficients  for  the  residuals  r2(t)  and  estimated 
coefficients  for  autoregressive  models  of  orders  p  =  1,2,3  and  4.  Note  that 
the  underlined  values  are  the  first  four  partial  correlations  for  the  data. 


Correlations 

I  Estimated  coefficients 

E. 

$(1)  =  0.8511 

aj  *  0.8511 

1 

p(2)  -  0.6973 

ax  =  0.9343; 

a„  =  -0.0975 

2 

p(3)  *  0.5718 

a^  a  0.9355; 

a2  =  0.1098; 

a3  »  0.0129 

3 

p(4)  *  0.4710 

=  0.955; 

a2  *  0.109; 

a3  =  0.007;  a*  =  0.0006 

4 

PERIODOGRAM  FOR  SECOND  RESIDUALS 
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ESTIMATED  AND  THEORETICAL  AUTOCORRELATIONS  FOR  SECOND  RESIDUALS 


FIGURE  13 


Consequently  the  first  residuals  are  pre-whitened  as  in  (4.7)  above,  using 
p  *  2  and  aj  a  0.9343  and  ^  s  -0.0975  derived  from  the  values  Ap(l)  ■  0.8511 
and  p(2)  *  0.6973.  The  resulting  third  residuals,  ra(t),  and  their  auto¬ 
correlations  are  shown  in  Figure  15  with  the  logarithm  of  the  normalized 
periodogram  in  Figure  16.  The  choice  of  frequencies  removed  on  the  basis  of 
the  first  residuals  is  confirmed.  To  explain  the  hump  in  the  autocorrelogram 
(Figure  14),  several  other  cyclic  components  were  removed,  e.g.  those  at 
p  =  19,  20  and  243.  These  changed  neither  the  autocorrelations  nor  the 
residual  variance  of  the  data. 

We  discuss  now  details  of  the  detrending  of  the  data.  Table  4  shows  the 
sanple  variance  of  the  original  data  (raw  data  minus  grand  mean),  the  first 
residuals  (removal  of  linear  trend),  the  second  stage  (removal  of  linear  trend 
and  s(t)  to  give  r2(t)),  and  the  fourth  residuals  which  will  be  discussed 
below. 


Stage 

Sample  Variance 

Table  4 

Sample  s.d. 

%  variance  removed 

Raw  data 

2.624 

1.62 

8.4% 

rx(t) 

2.403 

1.55 

51.9% 

r2(t) 

1.040 

1.02 

28.7% 

r4(t) 

0.288 

0.537 

10.9%  Remainder 

Note  that  about  half  of  the  variance  is  accounted  for  by  the  cyclic 
component  s(t).  A  break-down  for  the  0  =  13  cyclic  components  is  given  below 
In  Table  5.  Conponents  are  listed  in  descending  order  of  magnitude  of  their 
amplitudes,  .  The  variance  after  removal  of  each  of  these  (orthogonal) 
components  is  given  in  the  last  column.  The  annual  cycle  accounts  for  a  large 
part  of  the  variation,  namely  (. 81 )/l. 36  x  100%  =  59.6%  of  this  variance. 

These  estimated  coefficients  will  be  used  in  later  sections  in  a  simula¬ 


tion  study  of  the  model. 


THIRD  RESIDUALS  FROM  PRE-WHITENING[AR(2)] 
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Table  5 


j 

Tj (years) 

Yj 

8j 

Ao 

6j 

Variance 

2.40 

1 

12 

1.00 

-0.465 

-1.187 

1.275 

111.4° 

1.59 

2 

6 

2.00 

0.273 

-0.295 

0.402 

47.2° 

1.51 

3 

24 

0.50 

0.342 

-0.162 

0.378 

25.4° 

1.44 

4 

7 

1.71 

0.345 

0.112 

0.363 

342.0° 

1.37 

5 

3 

4.00 

-0.266 

0.211 

0.340 

218.4° 

1.32 

6 

9 

1.33 

-0.101 

0.302 

0.318 

251.6° 

1.28 

7 

2 

6.00 

0.074 

0.275 

0.285 

285.1° 

1.23 

8 

11 

1.09 

0.097 

0.250 

0.268 

291.2° 

1.19 

9 

18 

0.67 

-0.187 

-0.180 

0.260 

136.1° 

1.15 

10 

8 

1.50 

0.079 

-0.243 

0.256 

72.0° 

1.12 

11 

5 

2.40 

-0.176 

-0.174 

0.247 

135.3° 

1.08 

12 

94 

0.13 

-0.155 

-0.152 

0.218 

135.6° 

1.06 

13 

4 

3.00 

0.194 

-0.097 

0.217 

26.6° 

1.04 

6.  Examination  of  the  random  component,  e(t) 

The  second  residuals,  the  original  data  minus  the  linear  trend  and  the 
seasonal  trend  s(t)  containing  13  components,  are  whitened  (4.8)  under  the 
assumption  of  the  adequacy  of  an  AR (2)  model  to  give  the  fourth  residuals 


r4(t )  -  rx (t )  -  a^tt-l)  -  a2r2(t-2). 


The  first  use  of  these  residuals  is  as  a  test  of  the  adequacy  of  the 
modelling  which  has  been  done  up  to  this  point.  There  are  two  issues  which 
arise;  (1)  the  physical  adequacy  of  the  model, and  (2)  the  statistical  adequacy. 
Any  given  model  may  satisfy  either,  neither  or  both  of  these  requirements. 

For  physical  adequacy  we  require  that  the  model  describe  a  substantial  amount 
of  the  variation  in  the  data  with  as  few  terms  as  possible  (i.e, ,  the  principle 


of  parsimony),  and  hope  that  at  least  the  dominant  terms  have  oceanographic 
interpretations.  For  statistical  adequacy  we  have  made  the  assumption  that  the 
fourth  residuals  are  estimates  of  the  n(t)'s  and,  as  such,  should,  if  the  model 
is  reasonable,  be  approximately  uncorrelated. 

The  residuals  themselves  are  plotted  in  the  top  panel  of  Figure  17.  There 
is  an  anomaly  at  about  3100  days  associated  with  isolated  extremes  in  SSTs 
encountered  during  the  fall  of  1979.  It  is  in  fact  a  reminder  that  the  model  we 
have  derived  is  a  global  model  which  will  not  account  for  transients  or,  say, 
subtle  seasonal  effects.  Except  for  that  anomaly,  the  fourth  residuals 
generally  fall  in  the  range  of  ±  2.0  which  is  considerably  smaller  than  the 
range  of  the  original  data  (8.0  to  17.0).  Moreover,  in  Table  4  we  have  shown 
that  the  fourth  residuals,  which  represent  the  unpredictable  part  of  the  data, 
account  for  only  10.1%  of  the  total  variation.  Oceanographic  interpretations 
for  some  of  the  parameters  used  in  constructing  the  model  are  presented  later, 
as  well  as  simulation  results  which  show  how  well  the  model  mimics  in  behavior 
the  data. 

The  autocorrelation  function  for  the  fourth  residuals  is  given  in  the 
lower  panel  of  Figure  17.  The  bands  are  approximate  9b%  confidence  bands  for 
individual  correlations.  Again  there  may  be  an  anomaly  at  lag  20,  but  since 
this  is  a  maximum  of  100  estimated  correlations  it  probably  is  not 
significantly  large.  A  more  precise  test  that  the  r4(t)s  are  uncorrelated 
is  given  in  Figure  18. 

The  top  panel  is  the  unsmoothed,  normalized  periodogram  and  the 
bottom  panel  is  the  cumulative  periodogram.  Barely  distinguishable  upper 
and  lower  95%  and  99%  confidence  intervals  for  a  flat  spectrum  based  on 
Kolmogorov-Sml rnov  (KS )  statistics  are  shown.  Since  the  cumulative  periodogram 
lies  within  these  bands,  the  Durbin  test  for  a  flat  spectrum  (see  Jenkins  and 
Watts,  1969,  p.234-7,  or  Cox  and  Lewis,  1966,  p. 169 )  accepts  the  hypothesis  of 
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FIGURE  17 


F 


a  flat  spectrum.  The  value  of  the  KS  statistic  is  given  as  0.5136  and  compares 
to  the  upper  5%  point  of  the  statistic  which  has  value  1.355.  In  addition  the 
Anderson-Darling  test  statistic  (value  0.9139)  is  given;  this  test  statistic  is 
more  sensitive  to  departures  at  low  or  high  frequencies  (Cox  and  Lewis,  1966, 
p.169)  and  in  this  case  accepts  the  hypothesis  of  no  correlation. 

Thus,  we  have  established  the  adequacy  of  the  model.  Note  that  if  the 
e(t)  process  is  a  Normal  process,  then  a  flat  spectrum  implies  independence, 
but  not  otherwise.  For  this  reason,  and  in  order  to  complete  the  model,  we 
examine  the  distribution  of  the  fourth  residuals. 

Figure  19  shows,  in  the  top  panel,  a  crude  histogram  of  the  fourth 
residuals,  together  with  sample  statistics.  The  estimated  kurtosis  value  of 
3.61,  which  for  the  sanple  size  of  N  *  4380  is  significantly  different  from 
zero,  indicates  non-normality  in  the  residuals.  This  is  confirmed  in  the 
probability  plot  in  the  lower  panel.  This  is  a  plot  of  the  ordered 
residuals,  say  X^,  against  their  approximate  expected  value  *”l(1/N+l) 

where  *_1(.)  is  the  inverse  normal  probability  function.  The  curvature  in 
this  plot  and  the  outliers  indicate  that  the  distribution  of  the  e(t)’s  is 
somewhat  non-normal,  in  fact  pinched  in  the  middle  and  long  in  the  tails.  The 
distribution  appears  symmetric;  this  is  reflected  in  the  estimated  coefficient 
of  skewness  of  -0.00147. 

No  attempt  has  been  made  to  model  these  c(t)'s  in  more  detail.  Their 
standard  deviation  is,  from  the  last  line  of  Table  4,  estimated  to  be  0.537. 

The  autoregressive  parameters,  a:  and  a2,  for  the  AR(2)  process  are  estimated 
In  Table  3  to  be  0.9343  and  -0.0975  respectively  and  the  coefficients  of  s(t) 
are  given  In  Table  5.  The  long-term  trend  is  m(t)  *  10.9  +  0.000374t. 
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FIGURE  19 


In  Figure  8  we  have  given  a  local  smoothing  of  the  data  which  shows  the 
main,  gross  features  of  the  SSTs.  Although  there  are  many  ways  in  which  the 
data  can  be  smoothed  and  there  is  the  choice  of  an  appropriate  bandwidth,  the 
procedure  is  free  of  the  need  to  postulate  a  global  model.  In  fact  the 
smoothing  (Figure  8)  and  the  global  fit  of  s(t)  in  Figure  11  give  very 
similar  information.  Why  then  consider  global  models  of  m(t)  and  s(t)  and  the 
subsequent  modelling  of  e ( t ) ? 

The  answer  lies  in  our  objectives. 

( 1 )  We  want  to  use  the  model  to  i dent i fy  relevant  "scales"  in  the  data? 
The  question  of  scales  is  rather  nebulous,  with  different  answers 
coming  out  of  different  parts  of  the  model.  With  regard  to  s(t), 
there  is  clearly  a  "scale"  of  one  year  and  longer  scales  that  apply, 
for  example,  to  the  occurrences  of  El  Ninos.  Locally,  the  high 
correlation  in  e(t)  gives  another  scale,  measured  perhaps  by  the 
autocorrelation  function  of  the  e(t)'s  given  in  Figure  14.  And,  of 
course,  if  diurnal  effects  are  of  interest,  then  there  is  a  scale  of 
interest  shorter  than  the  one  day  sampling  interval  used  to  collect 
the  data. 

(2)  With  respect  to  scales,  we  want  to  know  how  often  should  the  data  be 
sampled?  This  is  a  question  which  cannot  be  answered  without  a 
specific  objective  in  mind.  For  instance,  if  we  want  information 
about  daily  variation,  we  need  to  sample  more  often  than  once  a  day. 
If  we  are  concerned  with  longer-term  variation,  we  can  relate  the 
sanpling  rate  to  the  accuracy  of  predictions.  One  component  of  the 
prediction  is  the  cyclic  and  long-term  trend;  how  accurately  are  they 


estimated  and  what  would  be  lost  by  sampling  once  every  two  or  three 
days?  In  Section  9  we  examine  sampling  variation  in  the  parameters 
in  m(t)  and  s(t)  by  generating  (simulating)  independent  realization 
of  SST  based  on  the  model. 

(3)  The  model  can  be  used  to  interpolate  missing  data.  This  is  not  a 
major  problem  in  the  Granite  Canyon  data  set  and  thus  it  will  not  be 
discussed.  However,  gaps  are  a  frequent  problem  with  many 
oceanographic  data  sets  (e.g.  in  the  Farallon  Island  data  mentioned 
previously). 

(4)  The  model  can  be  used  to  predict  or  forecast  future  temperatures, 
i.e.  given  the  data  up  to  today,  what  will  they  be  for  several 
subsequent  days.  We  show  in  the  next  section  how  the  model  is  used 
to  give  predictions  for  up  to  two  days  in  advance  of  current  data. 
These  predictions  will  be  shown  to  be  better  than  the  usual 
climatological  prediction  based  on  averages. 

If  we  attempt  to  use  a  climatological  model  (i.e.  an  average 
across  the  12  years  or  an  average  across  the  smoothed  values  for  12 
years)  to  address  the  questions  above,  we  sacrifice: 

(1)  The  variation  from  year-to-year,  so  that  events  like 

El  Ninos  and  the  long  term  trend  represented  by  m(t)  are 
masked  by  the  averaging. 

(2)  The  local  information  represented  by  e(t).  It 
seems  intuitively  reasonable  to  make  use  of  the  water 
temperature  today  If  we  need  to  make  a  prediction  for 
tomorrow.  Quantitatively,  Table  4  has  shown  that  this 
local  information  accounts  for  28.7%  of  the  variation  in  the 
data. 
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(3)  Decomposition  of  the  different  factors  affecting  sea-surface 
temperature,  retaining  only  a  picture  of  the  gross  behavior. 
Also  we  can  no  longer  consider  separately  the  yearly,  6 
monthly,  or  tidal  cycles. 

Hence,  we  find  the  model  to  be  more  useful  than  climatology  for 
summarizing  the  data.  We  use  the  next  two  sections  to  demonstrate  Its  value. 

8.  Forecast  of  SSTs  from  March  1,  1983  to  March  30,  1983 

We  now  forecast  SSTs,  one  and  two  days  ahead,  using  the  model  developed  in 
the  previous  sections.  Since  data  are  now  available  for  the  month  of  March, 
1983,  a  period  beyond  the  12  year  period  used  in  constructing  the  model, 
predictions  are  made  for  comparison  with  the  observations.  The  previous  twelve 
years  of  data  are  used  to  estimate  the  parameters  in  the  model.  These  forecast 
values  are  compared  to  the  actual  observations  and  to  a  climatological  forecast 
using  the  average  temperature  for  each  given  day  across  the  twelve  year 
ensemble. 

Given  an  autoregressive  process,  i.e,  the  second  residuals  r2(t)  as  in 
equation  (4.1)  of  Section  4,  with  known  coefficients  ax  and  a2, 

r2(t)  =  a^U-l)  +  a2r2(t-2)  +  e(t), 

and  that  the  i.i.d  residuals  are  normally  distributed  with  mean  zero,  then  the 
best  predictors  one  and  two  steps  ahead  are  obtained  by  setting  e(t+2)  and 
e(t+l)  equal  to  their  expected  values  of  zero.  We  obtain 


and 


Pred  (r2(t+l))  ■  a^t)  +  a2r2(t-l)  , 

Pred  (r2(t+2))  *  (a1z+a2)r2(t)  +  a^r^t-l)  , 


where  ^  and  a2  are  least  squares  estimates  of  al  and  a2  based  on  the 
observed  12  years  of  data.  For  periods  far  enough  ahead,  these  predicted 
values  would  eventually  converge  to  zero.  This  is  simpler  to  see  if  a2  =  0 

A  I 

so  that  we  have  an  AR ( 1 )  process.  Then  Pred(r2(t+k))  =  a*  r2(t)  . 

We  also  have  estimates  m(t)  and  s(t)  for  the  linear  and  seasonal  trends 
for  t=l  to  4380,  and  these  can  be  computed  beyond  that  range  to  give 
projections  m(t+l),  m(t+2),  and  s(t+l),  s(t+2).  The  forecasts  are  then 

Pred(y(t+1))  a  m(t+l)  +  s(t+l)  +  Pred(r2(t+1)) 

Pred(y(t+2))  =  m(t+2)  +  s(t+2)  +  Pred(r2(t+2)) 

These  forecasts  have  been  conputed  for  the  moving  twelve  year  period  out 
to  March  30,  1983,  together  with  the  climatological  forecast,  and  all  are  given 
in  Table  6. 
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Date 

Table  6 

Forecasts  of  one  month's  SSTs  at  Granite  Canyon 

One-step  Two-step 

t  Actual  y(t)  Forecast  Error  Forecast  Error 

One-step 

Climatol. 

Forecast 

Error 

March 

1, 

1983 

4381 

14.3 

14.0 

-0.3 

11.5 

-2.8 

4382 

14.4 

14.0 

-0.4 

13.7 

-0.7 

11.4 

-3.0 

4383 

14.2 

14.1 

-0.1 

13.7 

-0.5 

11.5 

-2.7 

4384 

14.5 

13.9 

+0.6 

13.8 

-0.7 

11.4 

-3.1 

4385 

14.3 

14.2 

-0.1 

13.6 

-0.7 

11.2 

-3.1 

4386 

14.4 

14.0 

-0.4 

13.8 

-0.6 

11.3 

-3.1 

4387 

14.2 

14.1 

-0.1 

13.7 

-0.5 

11.4 

-2.8 

March 

», 

1983 

4388 

14.8 

13.9 

-0.9 

13.8 

-1.0 

11.6 

-3.2 

4389 

14.9 

14.5 

-0.4 

13.6 

-1.3 

11.6 

-3.3 

4390 

14.7 

14.5 

-0.2 

14.1 

-0.8 

11.5 

-3.2 

4391 

14.8 

14.3 

-0.5 

14.1 

-0.7 

11.4 

-3.4 

4392 

14.0 

14.4 

+0.4 

14.0 

0.0 

11.6 

-2.4 

4393 

14.2 

13.7 

-0.5 

14.1 

-0.1 

11.1 

-3.1 

March 

14, 

1983 

4394 

14.5 

13.9 

-0.6 

13.5 

-1.0 

11.0 

-3.5 

4395 

13.6 

14.2 

+0.6 

13.7 

+0.1 

10.7 

-2.9 

4396 

13.5 

13.3 

-0.2 

13.9 

-0.4 

11.0 

-2.5 

4397 

13.8 

13.3 

-0.5 

13.2 

-0.6 

10.9 

-1.9 

4398 

13.7 

13.6 

-0.1 

13.2 

-0.5 

10.9 

-1.8 

4399 

13.7 

13.5 

-0.2 

13.4 

-0.3 

11.0 

-2.7 

4400 

13.7 

13.5 

-0.2 

13.3 

-0.4 

10.9 

-2.8 

March 

21, 

1983 

4401 

13.5 

13.5 

0.0 

13.4 

-0.1 

11.1 

-2.4 

4402 

13.8 

13.3 

-0.5 

13.4 

-0.4 

11.0 

-2.8 

4403 

13.2 

13.6 

+0.4 

13.2 

0.0 

11.0 

-2.2 

4404 

13.5 

13.0 

+0.5 

13.4 

-0.1 

11.0 

-2.5 

4405 

13.3 

13.4 

+0.1 

12.9 

+0.4 

10.8 

-2.5 

4406 

12.7 

13.2 

+0.5 

13.2 

+0.5 

10.6 

-2.1 

4407 

12.6 

12.6 

0.0 

13.0 

+0.4 

10.2 

-2.4 

March 

28, 

1983 

4408 

12.5 

12.6 

+0.1 

12.6 

+0.1 

10.2 

-2.3 

4409 

12.4 

12.5 

+0.1 

12.5 

+0.1 

10.2 

-2.2 

4410 

12.8 

12.4 

-0.4 

12.5 

-0.3 

10.5 

-2.3 

12.4  " 

Av 

-0.14 

Av 

-.38 

Av 

-2.7 

s. 

d  0.37 

s.d 

0.42 

s.d. 

0.45 

m.  s. 

e  0.40 

m.  s.e. 

,  0.57 

m.  s.e. 

2.74 

First  one  can  see  that  the  one-day-ahead  forecast  Pred(y(t+1))  and  the 
two-day-ahead  forecast,  Pred(y(t+2)),  are  usually  too  low  (i.e.  negatively 
biased),  the  average  differences  being  -0.14  and  -0.38  respectively.  The 
climatological  forecast  is  even  worse.  Of  course,  this  is  probably  due  to 
the  effect  of  the  current  El  Nino.  The  global  prediction  of  m(t)+s(t)  cannot 
quite  track  the  anomalous  warming,  although  the  autoregressive  residuals 
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compensate  to  some  extent.  The  climatological  forecast  is  quite  far  off  and  the 
comparison  because  of  the  El  Nino  is  probably  not  fair,  but,  nevertheless,  the 
model  given  in  this  paper  does  give  a  better  forecast.  Note  that  the  s.d.  of  the 
of  the  global  forecast  increases  as  one  goes  from  one  step  ahead  to  two  steps 
ahead,  as  does  the  bias.  This  is  to  be  expected  since  less  local  information  is 
available  for  the  two-day-ahead  prediction.  A  better  over-all  measure  of  the 
goodness  of  the  forecast  is  the  mean  square  error.  This  is  given  in  all  three 
cases  in  the  last  row  of  Table  6. 

9.  Simulation  results 

Here  we  use  the  model  for  the  Granite  Canyon  SSTs  to  generate  independent 
replications  of  the  time  series.  These  can  be  considered  as  members  of  an 
ensemble  from  such  a  process,  with  the  proviso  that  the  model  and  the  model 
parameters  used  to  generate  the  time  series  are  those  estimated  from  the 
original  Granite  Canyon  series. 

One  such  realization  is  shown  in  Figure  20.  The  top  panel  is  the  original 
series  and  the  bottom  panel  is  the  simulated  series.  Note  that  the  simulated 
series  is  generated  out  one  year  beyond  March  1,  1983  and  this  projection  is 
shown  in  Figure  21.  This  is  not  a  true  forecast,  except  insofar  as  the  trend 
m(t)  and  the  cycles  s(t)  account  for  a  major  part  of  the  variation  in  the 
Granite  Canyon  SSTs.  We  note  however,  that  the  one  year  simulation  (1983-84) 
does  contain  an  apparent  spring  transition  and  thus  introduces  a  note  of  reality 
into  the  results. 

Ten  independent  realizations  of  the  Granite  Canyon  SSTs  were  generated  and 
used  to  compute  estimated  parameter  values.  In  these  simulations,  the  "white 
noise"  sequence  n(t)  was  taken  to  be  i.I.d  Normal  random  variables,  despite  the 
departures  from  Normality  shown  in  Figure  19.  An  alternative  would  have  been  to 
"bootstrap"  the  process  by  using  a  random  permutation  of  the  r  (t)'s. 
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FIGURE  20 


In  Table  7  we  give  the  results  of  the  simulation  for  the  parameters  a  and  6 


in  the  long  term  trend,  for  the  serial  correlations  at  lags  one  and  two  in  the 
second  residuals  r2(t),  for  the  associated  coefficients  ax  and  a2  in  the  AR(2) 
model,  and  for  the  amplitudes  and  A2  for  the  two  dominant  frequencies  in 
s(t)  with  periods  of  one  year  and  two  years.  These  are  also  denoted  A  xand  A 
in  Table  5.  The  nominal  values  are  the  values  estimated  from  the  real  data  which 


were  used  to  generate  the  simulations. 


Parameter 

Table  7 

Simulation  results  -  sampling  every  point 

S.d. 

Nominal  value 

Min 

Max 

Average 

a 

10.9 

10.78 

10.95 

10.91 

0.045 

6 

0.000374 

0.000350 

0.000391 

0.000360 

0.000011 

P(l) 

0.8511 

0.8251 

0.8543 

0.8403 

0.0074 

P(2) 

0.6973 

0.6559 

0.7067 

0.6752 

0.0141 

0.9343 

0.890 

0.961 

0.929 

0.0187 

a2 

-0.0975 

-.136 

-0.078 

-0.106 

0.0191 

Ax  (lyr)  1.275 

1.217 

1.334 

1.278 

0.036 

A2  (2yr)  0.402 

0.360 

0.429 

0.401 

0.023 

Since  the  estimated  s.d.  of  the  slope  of  the  (linear)  long  term  trend 
m(t)  =  a+6t  is  0.000011,  it  is  clear  that  the  tilt  in  the  data  is  statistically 
significant.  Note  that  gross  approximations  give  the  s.d. 's  of  the  estimates 
of  p ( 1 )  and  p ( 2)  In  a  sample  of  size  4380  as  1/ (4380) =  0.015  which, 
given  statistical  fluctuations,  are  in  reasonable  agreement  with  the  estimated 
s.d. 's  of  the  serial  correlations  of  0.0074  and  0.0141.  Note,  however,  that  the 
estimates  of  p ( 1 )  and  p ( 2 )  have  averages  below  the  nominal  values.  This  reflects 
the  known  fact  that  estimates  of  serial  correlations  with  high  values  are  biased; 


Tables  7  and  8  are  virtually  identical.  For  independent  data  the  s.d. 's  in 
Table  8  would  be  approximately  /?  times  those  in  Table  7.  The  results  there¬ 
fore  indicate  that  day-to-day  correlation  in  the  data  is  so  high  (local  scales 
are  so  long)  that  sampling  every  day  does  not  give  any  Improvement  to  the 
sampling  of  long  scale  (one  year  and  up)  phenomena.  Such  phenomena  are 
effectively  summarized  by  the  model  through  s(t)  and  m(t),  exclusive  of  e(t). 
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In  Table  9,  results  are  given  for  estimation  with  sampling  every  third  day. 
The  standard  deviations  are  beginning  to  increase,  but  not  sufficiently  to 
justify  any  finer  sampling. 

10.  Oceanographic  aspects  of  the  data 

10.1  Spring  Transitions 

Abrupt  decreases  in  SST  have  been  observed  in  the  annual  temperature  cycle 
at  Granite  Canyon  during  the  spring  and  appear  to  signal  the  onset  of  upwelling 
along  the  central  California  coast.  These  dramatic  changes  in  temperature  occur 
in  at  least  6  out  of  the  12  years  considered  and  have  been  individually 
identified  and  tabulated  previously  in  Table  1.  From  the  table  it  is  seen  that 
these  transitions  typically  occur  in  March  and  last  for  about  a  week.  Also  the 
associated  temperatures  drop  from  about  13  to  about  IOC.  We  have  examined 
one  spring  transition  in  particular,  the  spring  transition  in  1980  starting  on 
March  11th,  which  is  shown  in  Figure  22.  During  this  transition,  temperatures 
dropped  about  4C  over  a  period  of  7  days.  As  indicated  previously  in  Figure  3, 
the  entire  water  column  experienced  this  transition  down  to  at  least  500  m. 
Alongshore,  currents  measured  at  the  same  location  reversed  from  poleward  to 
equatorward  in  direction  at  each  of  four  depths  between  113  and  510  m  at  the 
same  time  as  the  decreases  in  temperature  occurred.  The  synoptic-scale  winds  at 
this  location  however  did  not  become  upwelling  favorable  until  just  after  this 
event  occurred.  The  latter  observation  suggests  that  this  event  may  have  been 
caused  by  wind  forcing  at  another  location.  During  a  similar  spring  transition 
during  1981,  8rink,  Stuart,  and  Van  Leer  (1983)  also  concluded  that  local  wind 
forcing  was  not  responsible  for  the  transition  they  observed  off  Pt.  Conception 
during  that  year.  Finally  IR  satellite  imagery  acquired  a  few  days  before  and 
immediately  following  the  event  (not  shown)  indicated  a  change  from  isothermal 
surface  conditions  to  conditions  suggestive  of  intense  upwelling  in  the  vicinity 
of  Pt.  Sur. 
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ORIGINAL  GRANITE  CANYON  DATA- 


FIGURE  22 


The  relative  times  at  which  the  1980  spring  temperature  transition  arrived 
at  the  four  locations  along  the  California  coast  for  which  we  have  data  are  seen 
in  Figures  4,  5,  and  6.  This  event  is  seen  to  arrive  first  at  Diablo  Canyon, 
second  at  Granite  Canyon,  third  at  the  Farallons,  and  finally  at  Pacific  Grove. 
Phase  speeds  for  poleward  propagating  disturbances  along  the  coasts  of  North  and 
South  America  generally  range  between  about  35  and  280  Km/day.  In  particular, 
Allen  and  Romea  (1980)  indicate  that  phase  speeds  of  the  order  of  175  Km/day 
should  be  expected  for  barotropic  continental  shelf  waves  along  the  South 
American  coast  south  of  about  18°S.  We  take  this  value  as  representative  along 
the  North  American  coast  as  well,  since  the  continental  shelf-slope  regions 
along  both  South  and  North  America  are  generally  similar  (Shepard,  1973).  The 
difference  in  arrival  time  between  Diablo  Canyon  and  Granite  Canyon,  which  are 
about  175Km  apart,  is  approximately  a  day,  a  time  separation  consistent  with  the 
propagation  speed  of  175  Km/day  for  barotropic  continental  shelf  waves  in  this 
region.  However,  the  separation  in  arrival  times  between  Granite  Ca:<yon  and  the 
Farallons  is  considerably  greater  (~  6  days)  even  though  the  distance  between 
these  points  is  again  about  175  Km.  The  time  separation  in  this  case  is  not 
consistent  with  propagation  of  continental  shelf  waves  since  phase  speed  should 
in  fact  increase  due  to  the  slight  increase  in  the  Coriolis  parameter  and  an 
increase  in  width  of  the  continental  shelf  in  this  region.  Consequently  we 
conclude  that  these  results  are  not  necessarily  consistent  with  shelf  wave 
dynamics.  We  also  note  that  the  sea  surface  terrperatures  are  only  sampled  once 
a  day.  Thus  our  ability  to  resolve  the  phase  speeds  of  these  waves  is  limited. 

10.2  Annual  cycles 

A  mean  annual  temperature  cycle  has  been  calculated  for  the  Granite  Canyon 
data,  realizing  from  our  previous  statistical  calculations  that  a  mean  value  may 
not  be  a  truly  representative  measure  of  central  tendency  in  our  data.  Points 


corresponding  to  the  average  temperatures  and  a  smoothed  version  of  that  average 
(with  a  window  of  25  points)  are  both  presented  in  Figure  23.  According  to  List 
and  Koh  (1976),  and  Reid,  Roden  and  Wyllie  (1958),  a  well-defined  annual  cycle 
in  SST  might  not  be  expected  off  Granite  Canyon,  due  to  warming  during  the 
winter  from  the  Davidson  Current  and  cooling  during  the  spring  and  summer  from 
coastal  upwelling.  To  the  contrary.  Figure  23  shows  a  well-defined  annual  cycle 
at  Granite  Canyon,  with  a  range  of  about  3C.  However  the  annual  range  of 
variability  in  temperature  is  reduced  at  Granite  Canyon  compared  to  the  range  of 
variation  usually  found  at  this  latitude.  This  is  shown  in  the  insert  in  Figure 
23  where  deep  ocean  SSTs  (Robinson,  1976)  are  compared  to  the  Granite  Canyon 
SSTs.  The  effect  of  averaging  the  data  has  been  in  part  to  "smear"  the 
occurrence  of  the  spring  transition  over  a  two-month  period  from  about 
mid-February  to  mid-April.  Maximum  temperatures  are  seen  to  occur  during 
mid-October,  about  a  month  later  in  the  year  than  might  be  expected  at  this 
latitude.  Obviously  the  phase  relationships  between  such  factors  as  the  onset 
and  cessation  of  the  Davidson  Current,  and  the  onset  and  cessation  of  coastal 
upwelling  determine  in  part  the  detailed  structure  of  the  annual  temperature 
cycle  off  Pt.  Sur  in  any  given  year, 

Skogsberg  (1936)  originally  specified  three  oceanic  seasons  for  Monterey 
Bay,  an  upwelling  season  from  about  February  to  July,  an  oceanic  season  from 
about  July  to  November,  and  a  Davidson  Current  season  from  about  November  to 
February.  According  to  the  smoothed  average  annual  cycle  presented  in  Figure 
23,  at  least  two  "seasons"  can  be  identified,  an  upwelling  season  and  a 
non-upwelling  season.  Based  on  our  data  the  distinction  between  an  oceanic 
season  and  a  Davidson  Current  season  is  not  clear.  We  note  however  that 
Skogsberg  used  spatially  distributed  temperature  data  as  well  as  chemical  and 
biological  observations  to  construct  his  seasonal  description.  Also  his 


seasonal  model  was  derived  for  Monterey  Bay  and  as  a  result  it  may  not  be  valid 
to  extrapolate  his  model  to  regions  outside  the  Bay. 

10.3  El  Nino  events 

In  addition  to  the  seasonal  variation  in  temperature  at  Granite  Canyon, 
inter-annual  variability  is  also  evident,  as  mentioned  previously  in  Section  3. 
The  large  increases  in  temperature  seen  during  the  winters  of  1972,  1976,  and 
1983  correspond  to  major  El  Nino  warming  events  that  took  place  during  those 
years.  These  events  have  also  been  observed  off  central  California  through 
biological  measurements  and  measurements  of  sea  level  by  Chelton,  Bernal  and 
McGowan  (1982).  The  increase  in  temperature  seen  during  the  latter  part  of  1979 
is  not  associated  with  an  El  Nino  warming  event  (Enfield,  1983). 

Of  particular  interest  is  the  present  major  El  Nino  warming  event.  This 
warming  was  shown  to  extend  at  least  600  Km  off  the  coast  of  central  California 
during  February  1983  (Linder  and  Breaker,  1983).  The  effects  of  the  present 
warm  anomaly  are  clearly  evident  in  Figure  24  which  compares  the  3/1/82  to 
3/1/83  time  period  (i.e. ,  the  solid  curve)  to  the  12  year  average  data  (i.e. , 
the  dashed  curve).  Early  in  1982  warming  off  Granite  Canyon  appeared  to  be 
intermittent.  However,  by  late  1982  warming  at  this  location  had  become 
well-established.  As  of  February  1983,  the  warm  anomaly  at  Granite  Canyon  had 
reached  almost  3C.  A  smoothed  version  of  Figure  24  is  shown  in  Figure  2b  and 
helps  clarify  the  present  departure  in  SST  from  the  12-year  mean  at  this 
location. 

10.4  Interpretation  of  higher  frequency  variations 

From  an  oceanographic  viewpoint,  the  perlodogram  shown  previously  in 
Figure  9  is  somewhat  surprising  for  its  lack  of  any  major  spectral  components 
beyond  100  cycles.  For  example,  spectra  of  winds  along  the  California  coast 


often  show  higher  amplitudes  at  periods  between  3  and  10  days  (Van  Patten, 

1981).  There  are  three  possible  explanations  for  this  apparent  lack  of  con¬ 
sistency  between  wind  and  temperature  spectra.  First,  the  periodic  effects  in 
the  winds  may  be  real,  but  continual  stirring  of  the  near-surface  environment 
where  observations  of  ocean  tenperature  are  made,  may  inhibit  these  periodic  in¬ 
fluences  on  SSTs.  Alternatively,  higher  amplitudes  in  the  coastal  wind  spectra 
may  be  caused  by  correlation  in  the  wind  data  rather  than  by  underlying  period¬ 
icities.  As  such,  a  periodic-like  response  in  SST  to  the  local  winds  should  not 
necessarily  be  expected.  A  third  possibility  of  course  is  that  the  spectrum  of 
the  wind  field  at  this  particular  location  has  no  discrete  components  to  start 
with.  Data  to  check  this  third  possibility  unfortunately  do  not  exist  over  the 
full  12  year  period  but  there  is  evidence  from  other  locations  (Hugus,  1982) 
that  3  to  10  day  periodicities  do  not  show  up  consistently  in  wind  observations. 

Originally  the  possibility  of  aliasing  by  the  semiduirnal  and/or  diurnal 
tides  was  a  concern.  Analysis  of  similar  observations  at  other  locations  along 
the  California  coast  by  List  and  Koh  (1976)  indicated  that  tidal  aliasing  (or  a 
natural  periodicity  corresponding  to  the  spring/neap  tides)  by  the  semidiurnal 
tide  at  14.78  days  could  be  important.  To  the  contrary  we  do  not  find  evidence 
that  the  alias  of  either  tidal  component  (148  cycles  for  the  diurnal  tide  and 
296  cycles  for  the  semidiurnal  tide)  appears  in  the  calculated  periodogram. 

In  examining  the  periodogram  of  second  residuals,  particularly  at  high 
frequencies.  It  Is  of  interest  to  briefly  consider  the  likely  Influence  of 
turbulence  on  these  results.  Turbulence  production  undoubtedly  occurs  from 
the  breaking  of  incoming  surface  waves  and  perhaps  from  local  drift  currents  as 
well.  In  this  regard  a  portion  of  the  periodogram  of  second  residuals  has  been 
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smoothed  and  replotted  on  log  coordinates  in  Figure  26.  The  uniform  decrease 

in  spectral  amplitude  with  increasing  frequency  is  reminiscent  of  the  spectral 

decay  often  associated  with  the  inertial  subrange  from  homogeneous  turbulence 

theory  (Batchelor,  1963).  Between  21)0  (~  2  days)  and  2U00  cycles  (~22  days), 

the  spectral  slope  is  essentially  constant  at  -9/6,  close  to  the  value  of  -6/3 

predicted  by  classical  turbulence  theory  (Phillips,  1977).  According  to  this 

theory  temperature  is  assumed  to  be  a  conservative  quantity,  an  assumption  which 

could  be  quite  unrealistic  at  the  sea-surface.  In  our  case  time  replaces 

distance  which  is  normally  taken  as  the  independent  variable,  and  if  it  is 

assumed  that  these  coordinate  frames  are  linearly  related,  then  the  spectral 

slopes  observed  here  are  directly  comparable  with  those  usually  derived  in 

oceanic  turbulence  calculations  as  a  function  of  distance.  In  turbulence  models 
* 

where  temperature  is  not  assumed  to  be  conservative,  spectral  slopes 
considerably  greater  than  -6/3  are  predicted  (Panchev,  1980).  Measurements  of 
the  spatial  surface  temperature  spectrum  by  Saunders  (1972)  and  Holladay  and 
O'Brien  (1976)  for  example,  yielded  spectral  slopes  of  -2.2  and  -3  respectively, 
suggesting  that  temperature  should  not  be  treated  as  a  conservative  quantity  at 
least  under  certain  conditions. 

In  conclusion,  our  earlier  analyses  suggest  that  the  (likely)  influence 
of  turbulence  in  our  data  corresponds  to  the  autoregressive  component  in  our 
statistical  decomposition  since  the  spectrum  of  the  fourth  residuals  is 
essentially  flat.  A  more  detailed  analysis  of  our  data  would  include  a 
comparison  of  the  spectra  for  various  orders  of  autoregression  with  spectra 
predicted  from  various  theoretical  models  of  ocean  turbulence. 
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